knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE)
knitr::knit_meta(class=NULL, clean = TRUE)
## list()
if (!require("FactoMineR")) install.packages("FactoMineR", dependencies = TRUE) ; library(FactoMineR)
if (!require("factoextra")) install.packages("factoextra", dependencies = TRUE) ; library(factoextra)
if (!require("plotly")) install.packages("plotly", dependencies = TRUE) ; library(plotly)
if (!require("tibble")) install.packages("tibble", dependencies = TRUE) ; library(tibble)
if (!require("gridExtra")) install.packages("gridExtra", dependencies = TRUE) ; library(gridExtra)
#install.packages('RMySQL', repos='http://cran.us.r-project.org')
#if (!require("")) install.packages("")
if (!require("backports")) install.packages("backports", dependencies = TRUE) ;
if (!require("curl")) install.packages("curl", dependencies = TRUE) ;
if (!require("haven")) install.packages("haven", dependencies = TRUE) ;
if (!require("openxlsx")) install.packages("openxlsx", dependencies = TRUE) ;
Charge the libraries for this analysis
read files, with column names already in the file (headers=T), and data separated by colon (sep=";")
bacteria.df1 <- read.csv("./2020092301-MJ-raw-Bacteria.csv", header=T, row.names=NULL, sep=",")
maxima1.df1 <- read.csv("./2020092301-MJ-raw-Maxima1.csv", header=T, row.names=NULL, sep=",")
Column names of bacteria.df dataframe
colnames(bacteria.df1)
colnames(maxima1.df1)
Select the columns to be kept
bacteria.newColumns <- c("NAME.id", "EXPERIMENT.name", "MAXIMA.Maxima1.count",
"SHAPE.angularity.mean", "SHAPE.area", "SHAPE.circularity",
"SHAPE.curvature", "SHAPE.length", "SHAPE.width",
"SHAPE.perimeter")
maxima1.newColumns <- c("NAME.id", "EXPERIMENT.name", "INTENSITY.ch1.mean", "INTENSITY.ch2.mean",
"INTENSITY.ch3.mean",
"PARENT.location.cartesian.x", "PARENT.location.cartesian.y",
"PARENT.location.orthogonal.x", "PARENT.location.orthogonal.y")
bacteria.df2 <- bacteria.df1[,which(colnames(bacteria.df1) %in% bacteria.newColumns)]
maxima1.df2 <- maxima1.df1[,which(colnames(maxima1.df1) %in% maxima1.newColumns)]
bacteria.df3 <- bacteria.df2 %>% separate(EXPERIMENT.name, c("DATE","Time","Order","Phage"))
maxima1.df3 <- maxima1.df2 %>% separate(EXPERIMENT.name, c("DATE","Time","Order","Phage"))
Visualize the distribution of foci across the cells at different times
max.df <- maxima1.df3
InfTime <- c("T0","T10","T20","T30")
PhagStrain <- c("WT", "D257", "DR", "DLDR")
maxPlot <- function(i,j,k){
test <- max.df[(max.df$Phage==i) & (max.df$Time==j),]
x <- test$PARENT.location.orthogonal.x
y <- test$PARENT.location.orthogonal.y
ggplot(test, aes(y, x))+geom_hex(bins = 25)+scale_fill_gradientn("", colours = rev(rainbow(10, end = 4/6)))+
scale_y_continuous(breaks=seq(-5,5,2), limits=c(-5,5))+
scale_x_continuous(breaks=seq(-2,2,1), limits=c(-2,2))+
annotate("text", x=0, y=4.5, label=k, size=2.5) +
theme_classic()+
theme(legend.position = "none",
axis.text.x = element_blank(), axis.title.x = element_blank(),
axis.text.y = element_blank(), axis.title.y = element_blank())+
geom_ellipse(aes(x0 = 0, y0 = 0, a = 1, b = 3.5, angle = 0, m1=2, m2=3))
}
p <- list()
for(i in PhagStrain){
for (j in InfTime){
k <- paste(i,j,sep=" - ")
p[[k]] <- maxPlot(i,j,k)
}}
g <- do.call(grid.arrange, c(p,ncol=8))
#ggsave("./Maxima1-orthogonal.png", g, width = 16, height = 15, units = "cm" )
#Distribution of Parameters
colnames(bacteria.df3)
## [1] "NAME.id" "DATE" "Time"
## [4] "Order" "Phage" "MAXIMA.Maxima1.count"
## [7] "SHAPE.angularity.mean" "SHAPE.area" "SHAPE.circularity"
## [10] "SHAPE.curvature" "SHAPE.length" "SHAPE.perimeter"
## [13] "SHAPE.width"
bacteria.df3$Phage <- factor(bacteria.df3$Phage,levels=PhagStrain)
maxima1.df3$Phage <- factor(maxima1.df3$Phage,levels=PhagStrain)
mycolors = c(brewer.pal(name="Dark2", n = 8), brewer.pal(name="Paired", n = 6))
Distribution of area
plot_ly(bacteria.df3, x=~Phage, y=~SHAPE.area, color=~Time, colors = mycolors, type="box",
text = paste(bacteria.df3[,"NAME.id"])) %>% layout(boxmode = "group", xaxis = list(title = "Phage"), autosize = F, width = 750, height = 500, automargin = TRUE)
Plot the relative DAPI intensity across time
font1 <- list(
family = "Arial, sans-serif",
size = 12,
color = "black"
)
plot_ly(maxima1.df3, x=~Phage, y=~INTENSITY.ch1.mean, color=~Time, colors = mycolors,
type="box", text = paste(maxima1.df3[,"NAME.id"])) %>%
layout(boxmode = "group",
xaxis = list(title = "Phage", titlefont=font1, tickfont=font1),
yaxis=list(title = "DAPI Intensity", titlefont=font1, tickfont=font1),
legend=list(font=font1),
autosize = F, width = 1000, height = 500, automargin = TRUE)
Build grid of p-values from pairwise comparison: 1. Function to perform Tukey HSD test:
TukeyAllParams <- function(paramToTest, df){
df.lm <- lm(df[,paramToTest] ~ Phage , data = df)
df.av <- aov(df.lm)
tukey.test <- TukeyHSD(df.av, conf.level = 0.9999)
params.df <- data.frame(Comparisons=rownames(tukey.test$Phage), tukey.test$Phage[,"p adj"])
params.df2 <- params.df %>% separate(Comparisons, c("Gene1", "Gene2"))
params.df3 <- params.df2[which((params.df2$Gene1 %in% RefPhage) | (params.df2$Gene2 %in% RefPhage)),]
params.df3[params.df3 == RefPhage[1]] <- NA
params.df4 <- data.frame(t(apply(params.df3,1,function(x) x[!is.na(x)])))
colnames(params.df4) <- c("Phage", paramToTest)
rownames(params.df4) <- NULL
return(params.df4)
}
RefPhage <- c("WT")
bacteria.newCol2 <- colnames(bacteria.df3)[6:length(colnames(bacteria.df3))]
PhageList2 <- setdiff(levels(bacteria.df3$Phage), RefPhage)
p.All <- data.frame(Phage=PhageList2)
for (i in bacteria.newCol2){
p.All <- merge(p.All, TukeyAllParams(i, bacteria.df3[bacteria.df3$Time=="T10",]), by = "Phage", all.x = TRUE)
}
p.All2 <- column_to_rownames(p.All, "Phage")
p.All3 <- p.All2 %>%
rownames_to_column() %>%
gather(colname, value, -rowname)
p.All3$rowname <-factor(p.All3$rowname, levels=PhageList2)
p.All3$value <- cut(as.numeric(p.All3$value),breaks =c(-Inf, 1e-6, 1e-3, 0.01, 0.05, Inf))
p.plot <- ggplot(p.All3, aes(x=rowname, y=colname, fill=value))+geom_tile()+
theme(axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(color = "black", size=12),
axis.text.y = element_text(color = "black", size=12),
legend.text = element_text(color = "black", size=12))+
scale_fill_brewer(palette = "OrRd", direction=-1)
p.plot
#ggsave("./pVal-TukeyHSD-T30.png", p.plot, width = 20, height = 10, units = "cm" )